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Abstract 

A single, augmented Extended Kalman Filter (EKF), which simultaneously and autonomously 
estimates spacecraft attitude and orbit has been developed and successfully tested with real 
magnetometer and gyro data only. Because the earth magnetic field is a function of time and 
position, and because time is known quite precisely, the differences between the computed and 
measured magnetic field components, as measured by the magnetometers throughout the entire 
spacecraft orbit, are a function of both orbit and attitude errors. Thus, conceivably these 
differences could be used to estimate both orbit and attitude; an observability study validated this 
assumption. The results of testing the EKF with actual magnetometer and gyro data, from four 
satellites supported by the NASA Goddard Space Flight Center (GSFC) Guidance, Navigation, 
and Control Center, are presented and evaluated. They confirm the assumption that a single EKF 
can estimate both attitude and orbit when using gyros and magnetometers only. 


I. INTRODUCTION 

Many future low-budget missions, such as the NASA Small and Mid-Size Explorer 
Series and university class explorers, are looking for inexpensive and autonomous approaches to 
orbit and attitude estimation. One sensor being considered as a prime is the magnetometer 
because it is reliable and low cost. For these reasons, magnetometers have been the focus of 
several recent studies. Emphasis has been placed on using only the magnetometer to separately 
estimate the spacecraft trajectory 12 ' 3 and attitude 4 5 . Studies have also been performed which 
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show remarkable accuracy of the magnetometer in estimating attitude when accurate rate 
information is available' 1 . 

In this work, we present the design and test results of a single extended Kalman filter and its 
application to actual data obtained from four satellites; namely, Compton Gamma Ray 
Observatory (CGRO, NORAD No. 20580), Rossi X-Ray Timing Explorer (RXTE, NORAD No. 
23757), the Earth Radiation Budget Satellite (ERBS, NORAD No. 15354), and the Total Ozone 
Mapping Spectrometer-Explorer Platform (TOMS-EP, NORAD No. 23940). The present work is 
an extension of the research reported on in Ref. 7 where both attitude and trajectory have been 
successfully estimated using simulated magnetometer and rate data. In that work large initial 
errors were applied with a resulting accuracy of 4 km root sum square (RSS) in position and less 
than 1 degree RSS in attitude. 

In the work presented here, we summarize the derivation of the EKF, highlighting the 
development of the measurement matrix, which constitutes the crucial element in combining the 
dependence of the magnetic field residuals on both attitude and trajectory estimation errors. The 
satellites considered here vary in inclination (23 degrees for CGRO to 97 degrees for TOMS-EP) 
and in size of the magnetometer quantization error (0.3 to 6.4 milliGauss). Comparisons are 
made with operational estimates of position, velocity, and attitude that were very accurate with 
respect to the filter generated estimation errors (e.g. CGRO position measurement error was less 
than 1000 meters and that of ERBS was within 200 meters, etc.). These reference values were 
computed by NASA Goddard Space Flight Center (GSFC) flight dynamics personnel. We show 
the ability of the filter to overcome large initial errors and a comparison of the final accuracy 
achieved with each satellite. This work could prove valuable as a prime trajectory and attitude 
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estimation system for satellites with coarse accuracy requirements, or as a backup to a prime 
system in satellites with greater accuracy needs. 

II. EXTENDED KALMAN FILTER (EKF) ALGORITHM 

For clarity of the developments that will follow we review the special features of the EKF 
algorithm that was used in this work. The EKF algorithm is based on the following assumed 

models: 


System Model: 

X(t) = f(X(t),t) + w(t) (la) 


Measurement Model: 


y k+1 -hk+i(X( t k+i)) + v k+i 


(l.b) 


where w(t) is a zero mean white process noise, v k+1 is a zero mean white sequence measurement 
error, and X(t) is the state vector. In this work X(t) is defined as 


X T =[a, e, i, Q, to, 6 , C d ,b T , q T ] 


( 2 ) 


The first six elements of X(t) are the classical Keplerian orbital elements 8 ’ 9 which determine the 
spacecraft orbit, position and velocity; namely, the semi-major axis (a), eccentricity (e), 
inclination (i), right ascension of the ascending node (Cl), argument of perigee (©), and true 
anomaly (0). C d is the drag coefficient, b is a vector of the on board gyro constant drift rates, 

and q represents the attitude quaternion. 

The EKF used in this work is a unified filter that consists of an orbital part and an attitude 
part. Thus Eqs. (1) are written as 


X' 


"fo(X 0 (t),t)' 

4 - 

X" 

X 


_f,(X a (t),t)_ 


_ W a_ 


(3. a) 
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and 

=h k . 1 (X, ) (t kM ),X J (t kM )) + v k , I (3 b) 

where 

X 0 r =[a,e,i,Q,a),0,C a ] (4.a) 

and 

Xj=[b T ,q T ] (4 b) 

The processes, which are applied in order to derive the linearized dynamics and measurement 
models used in the EKF, treat the orbital and the attitude parts differently according to the 
specific physical natures of each part. The difference in treatment is expressed in the 
development of the measurement model as well as in the development of the dynamics model. 
The development of the orbital part of the measurement model is done by the customary partial 
differentiation, 3h/dX 0 , which yields a Jacobian matrix 10 whereas, as will be explained later, the 

attitude part is done using the perturbation method. Similarly, in the development of the orbital 
part of the dynamics matrix, the customary partial differentiation yields the Jacobian matrix 
df 0 /dX 0 . However, as will be shown later, the development of the dynamics model that 
describes the evolution of the attitude error between measurement updates is done entirely 
different. 

The state vector usually consists of three kinds of states; those that describe the dynamics 
of the system, those that are augmented due to the nonwhite nature of the stochastic processes 
that drive the dynamics model, and those that are augmented due to the same in the measurement 
noise. In our case the states a, e, i, Q, co, and 0 that appear in X of Eq. (2) are states that describe 

the system dynamics. The state C d is a non-white state driving the dynamics model, which was 
augmented via the state augmentation procedure. Similarly q is a state vector that describes the 
dynamics whereas b, the gyro bias state vector that drives q, is non-white, and as such is also 
augmented into the state vector. Sometimes it is also necessary to describe the measurement error 
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as a combination of white noise and colored noise 1 , which, as mentioned before, results in further 
augmentation of the non-white vector into the state vector. However in our filter this was not 
necessary. 

We now proceed to describe the measurement update stage and the propagation stage of 
the EKF algorithm used in this work. 

II.l EKF Measurement Update 

The measurement update of the state estimate and of the estimation error covariance are 
performed, respectively, as 

X k+ ,(+)*X k+1 (-)+K k+1 [y k+1 -h k+1 (X k+I (-))] (5.a) 

P k+1 (+) = [I- K k+1 H k+1 ]P k+1 (-)[I - K k+1 H k+1 ] T + K. k+1 R k+1 K k+1 (5.b) 

where (-) denotes the a-priori value and (+) denotes the a-posteriori value, and K k+I is the 
Kalman gain computed according to 

k w = wX,[h m p w (-a + R M r' <«> 

where 

H M = [HJH,1 1>1 (7) 

and where H 0 is the orbit part of the measurement matrix, and H a is the attitude part. The matrix 
P k+1 is the estimation error covariance matrix, and R k+1 is the covariance matrix of the zero mean 
white sequence v k+1 

Let us denote y k+I - h k+l (X k+l (-)) of Eq. (5.a) by z k+l . The latter is the data, which is used 
by the EKF to update the a-priori state estimate. This data, which is called effective measurement. 
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is the difference between the measurement and the estimate of the measurement. Therefore in our 
case the data, z k+l , used by the filter is computed as 


Z k-1 — ®(^k*l( 


( 8 ) 


here B mk+1 is the magnetic field vector measured by the magnetometer, B(X k+1 (-),t k+1 ) is the 

estimated magnetic field vector as a function of the estimated state at time t k+ , computed using a 

10 th order International Geomagnetic Reference Field (IGRF) model. 

Define x k+1 as the difference between the true state vector X k+1 and its estimate; that is, 

= [5a, 5e, 5i, SQ, 5co, 50, 5C d , 5b T , o T ] k+1 (9) 

A comparison between X, defined in Eq. (2), and x, defined in the last equation, reveals that 
other than q, the latter is a perturbation of the former. In x the expression for q is replaced by a . 
This is done because, as will be shown later, in the EKF update stage we express the attitude 
error by a vector of small angles, a , whereas in the propagation stage we express the attitude by 
the quaternion q. We use z k+1 , [Eq. (8)], to update the estimate, x k+1 of x k+1 as follows 

*k+l = Kk+l Z k+l (10) 

The updated vector, x k+1 , is used to update the a-priori state estimate, X k+1 (-) , as follows 

i,.,) 0 0 

The nature of g will be explained shortly. 


Development of the measurement matrix 

Figure l shows the geometry of the magnetic field vector. The U coordinates are earth 
fixed coordinates and the I are inertial coordinates. Both coordinate systems are fixed to the earth 
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center. The Magnetic field vector is expressed in the magnetic spherical coordinates, F, as 
follows B h r = [B r , B, m , B^jj], Begin the development of the measurement matrix using the 
perturbation technique, by writing the estimated magnetic field at the spacecraft location as 

B(X,t) = D'D^B F + v' (12) 

where B F is computed using the IGRF magnetic field model, the estimated position and time, 

A P 

D I is the estimated transformation matrix that transforms from magnetic spherical to inertial 

coordinates, D|, is the estimated attitude matrix which transforms the estimated magnetic field 
vector from the inertial to the body coordinate system, and v’ is the magnetic field model error. 
(Note that we consider the situation at time point t k+1 ; however we drop the subscript k+1 for 
clarity). On the other hand, the measured magnetic field vector, as measured by the 
magnetometer, can be written as 

B m = DbD, F B F + v m (13) 

where B F is the correct magnetic field vector in magnetic spherical coordinates, Df is the 
correct transformation from the magnetic spherical to the inertial coordinates, D{, is the true 
transformation from inertial to body coordinates, and v m is the magnetometer measurement 
error. 
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The effective measurement, z, is defined as follows. 


z = B m - B(X(-), t) = D{,D[B f + v m - D^DfBp - v' 


Write the transformation of B F to the estimated body coordinates as 

D{,DfB F = D[ ) DfB F + A(D{ ) DfB F ) (15) 

where A(DbDfB F ) is the error in the estimation of the magnetic field vector resolved in the 
body axes. Define 

v = v m - v' (16) 


Using the last two equations, Eq. (15) becomes 

z = A(D'D'B f ) + v (17) 

Now 

A(D‘Dj : B F ) = AD{,(DfB F ) + D{,A(D| : B F ) (18) 

The second term on the right hand side of Eq. (18) is the part of the effective measurement 
caused by errors in the orbital states. This term can be developed into H 0 x 0 where H 0 is the 
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measurement matrix for the orbital states and x 0 denotes the orbital error states, which are the 
first seven elements of x, that are defined in Eq. (9). The derivation of H o x 0 is given in Ref. 10 
where H 0 was developed using partial differentiation 


„ dh(X„,X ) 

H 0 = ° a . (19) 

dX 0 X = X 

The expansion of the first term leads to the measurement matrix for the attitude states. Rewrite 


that term as 


AD'(D, f B f ) = AD'B i 


( 20 ) 


where B, is the computed magnetic field vector in inertial coordinates. We define AD|, , the 

error in the transformation, as the difference between the correct body coordinates and the 
computed (estimated) body coordinate system. The computed coordinates are determined by the 

computed (estimated) transformation matrix, D{, , which is the transformation from the inertial to 
the computed (rather than the correct) body coordinate system. We write it as 

61 = dJ = dJdJ, (2i) 

SO 

ADt-Dt-D^DX-Dj, (22) 

For small attitude errors we can assume that the matrix is composed of small angles, thus 

"0 — \j/ S' 

D£=I- vj/ 0 — cp = I-[ax] (23) 

-9 (p 0 

where [ax] is the cross product matrix of a defined in the last equation, and where 
a T = [<p, 9, ip] . The vector a is a vector of three small Euler angles describing the attitude 
difference between the true and estimated attitudes. These errors are defined about the body x, y, 
and z axes respectively. Using Eq. (23) in Eq. (22) yields 
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AD;, = ([-[ax])Di- D b = -[a x ]D b 


(24) 


Substituting Eq. (24) into the first term on the right-hand side of Eq. (18) yields 

ADUDfB F ) = 4ax]D‘B, . -[ax]B b =[B b x]a (25) 

Following the last equation we substitute the first term of Eq. (18) by [B b x]a and, following 
earlier discussion, the second term by H 0 x 0 . Substituting the result into Eq. (17) yields 

z = [B b x]a + H 0 x 0 +v (26) 

The error state, x, is composed of the orbital error states denoted by x 0 , the error in estimating 
the gyro bias, 8b , and the vector of small angular errors in the attitude, a (see Eq. 9). Therefore, 
in order to relate z to x we need to introduce a matrix that relates 8b to z. However, 8b has no 
direct influence on z, therefore we can write Eq. (26) as 

z=[H o |0 3x3 |[B b x]]x + v (27) 

where 0 3x3 is a 3x3 zero matrix. (Note that b is introduced into the state vector, X, because it is 

a non-white vector driving the attitude. In order to comply with the requirement of the KF that 
the driving force be white, b is augmented with the original state vector. This is a well-known 
standard procedure. The existence of b in X results in the existence of 8b in x.) Following Eq. 
(27), the combined measurement matrix is given as 

H = [H o |0 3x3 |[B b x]] = [HJHJ (28) 

where H a = [0 3x3b | [B„x]]. Because B b is not known, B m , the magnetic field vector measured 
by the magnetometer, is used instead, thus 

H,=[0 M |[B„xj] (29) 

An explicit expression for H 0 is given in Ref. 10 using partial differentiation. 
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We realize that a new state vector, a, has emerged in the development of the 
measurement matrix. This vector of angular errors expresses the attitude error between the true 
and estimated attitude. W r e will show shortly that this vector will easily blend with the orbital 
error states. Using Eqs. (26) and (27), the effective measurement at time t k+1 is written as 

^k+l — ^k+l x k+l ^k+I (30) 

where 

= [8a, 5e, 5i, SQ, 5co, 50, 5C d , 5b\a T ] k+I (31) 

Define 

x L*i = [ 5a > 5e » Si, 8Q, 5©, 50, 5C d ] k+1 (32) 

from Eqs. (31) and (32) it is obvious that 

x Li = [*I,k+i>8b T , a k+1 ] (33) 

Similar to x d>k+1 , define 

Xl, k+ i = [a, e, i, Q, to, 0, C d ] k+1 (34) 


then from Eqs. (2) and (4) we obtain 

XL=[X 0 T k+ ..b T , q k+1 ] 


(35) 


With Eq. (35) we can now explain the nature of g in Eq. (11). Indeed, the update of the state 
vector estimate is a function of both X k+1 (-) and x k+I which is computed using Eq. (9). While 
the orbital and the gyro bias parts of X k+1 (-) are updated additively as follows 


Xm(+)" 

= 

XwH 

+ 

X o,k+l 

. b M (+) . 


. b w H . 


. 5 ^k + i. 


(36) 


the attitude part is updated multiplicatively as described next. It is well known that the three 
small angles in the vector, a, can be used to construct a quaternion dq as follows 11 

dq T =[!«£, 10, lvj>, l] = [i« T ,l] (37) 
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where dq is the estimate of dq which expresses the small attitude difference between our best 
estimate of the attitude and the correct attitude. In other word, if one rotates the correct 
coordinate system, whose orientation is expressed by q k ., , through the rotation expressed by dq , 
one reaches the current (a-priori) estimate of the orientation expressed by q ktl (-)- This is 
represented by 

qk + i(-) = qk + i®dq k+1 (38) 

where 0 denotes a quaternion product. Post-multiplying Eq.(38) by the inverse of dq k+l yields 

qk + i=4k + t(-)® d <Ik!i (39.a) 

If dq k+1 is known precisely then the operation q k+l (-) 0dq k J., yields q k+1 , but because we have 
only an estimate, dq k+l , of dq k+1 the product q k+1 (-) 0 dq k ^, yields only the estimate of q k+1 . 
Because dq k+ , is the a-posteriori estimate of dq k+1 , this new estimate of q k+1 is the a-posteriori 
estimate. In other words 


qk + i(+) = qk + i(-)® d qki. 


(39.b) 


The last equation can be viewed as a correction to q k+1 (-) which yields q k+1 (+). In order to 
assure the normality of the quaternion, normalization is performed at this stage, which yields 

4,.,(+) = 4,.,(+)1q,.,(+)| (39.c) 

Following Eq. (35) the full updated state vector is given by 

XL(+) = [X:. k .,(+),bL,qL(+)] (40) 


II. 2 EKF Propagation 

Recall the dynamics equation (3. a) 


-x/ 


■f 0 (X 0 (t),t)' 

4* 

* W o" 



.^»(Xj(t), t)_ 


- w a- 


(3. a) 
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The propagation of the covariance is given as 


P k+ 1 (-)=A k (X k ( + ))P k (+)A k r (X k (+)) + Q k (41) 

where Q k is the spectral density matrix of w(t) and A k is the approximated transition matrix. 
A k is computed using the following first order Taylor series expansion 

A k = I + F k • AT (42) 

where AT is the time interval between gyro measurements. As described earlier, the dynamics 
matrix of the orbit part is obtained by computing the Jacobian F 0 k as follows 


F 0 ,k = 


ggXpj t),t) 

ax 


X = X k (-) 


(43) 


The actual computation of F 0 k is based on the orbital dynamics equations, given in the ensuing. 

Based on Eq. (1), the propagation of the state estimate, is performed by solving the differential 
equation 

i(t) = f(X(t),t) (44) 

The updated estimate of the state vector, X k (+) serves as the initial condition of the solution, 

and the final solution at time t k+ i is X k+1 (-). The dynamics equation (44) consists of the orbital 
dynamics equation, which in our case, is non-linear, and the attitude kinematics, which is linear. 
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II. 2.1 Nominal Dynamics 


Orbital Dynamics 

The orbital dynamics describe the motion of a mass point in a central force field 
including drag 1 '. There are several possibilities to describe the orbital dynamics. The most 
natural possibility is by position, velocity and acceleration expressed in a Cartesian coordinate 
system. The advantage of this choice is in its ability to handle orbits which are completely 
circular but its disadvantage is in the fast change of the state that describes the dynamics. This in 
turn requires an increased computational load. One can also choose equinoctial variables 13 to 
describe the orbital dynamics. This choice too enables the handling of completely circular orbits. 
Moreover, all but one variable of the equinoctial variables change slowly [see Ref. 9, p. 143]. 
Finally, one can use Keplerian parameters, which, similarly to the equinoctial variables, include 
only one fast varying parameter. Although these parameters cannot describe perfectly circular 
orbits, as will be seen later, they can handle nearly circular orbits. For this reason and because 
Keplerian parameters are widely used and well understood we chose them for the present study. 

The differential equations for each of the orbital elements of the state vector X, including 
the drag coefficient are given as 10 

2a f* 

x,=a = y(2a-r)^ (45. a) 


2f. f r 

x, =e = — (cos0 + e) + — -sin0 
V V a 


: rf h . 

X, = 1 = — — COS0 


x, = Q = 


rf h 

h sini 


sin0 


(45. b) 


(45. c) 


(45. d) 
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(45.e) 

(45.0 


x 7 = C d = 0 (45. g) 

where f t , f n and f h are the along track, radial, and cross track perturbing accelerations due to the 
effect of drag, and where 


e‘ = 

CO +0 

(46. a) 

h = y 

W a (l~e 2 ) 

(46.b) 

a(l-e 2 ) 


r = 

l+ecos0 

(46. c) 


i, ,i i , 

(46. d) 


and n E is the earth gravitation constant. 


Attitude Dynamics 

The linear differential equation governing the dynamics of the attitude quaternion is given by 

q = ^q (47. a) 

where Q is a 4x4 skew-symmetric matrix containing the elements of the spacecraft rate vector 
as measured by the gyros 11 . We assume that the gyros have bias, and using the state 
augmentation approach 14 , mentioned in Section II, the three bias states were included in the state 
vector (Eqs. (2) and (4.b)). The bias states model is 

b = 0 (47. b) 
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It should be mentioned here that although the correct modeling of bias is as shown in Eq. (47. b), 
we added white driving noise for filter stability'^. 

II.2.2 Error-State Dynamics 
Orbital Error-State Dynamics 

The error dynamics model can be written as follows 



(48) 


where F 0 is the orbital part and F a is the attitude part of the error state dynamics. From Eq. (32) 
it is obvious that the dimension of F 0 is 7x7. The elements of this matrix are obtained when 
applying the partial differentiation of Eq. (43) to Eqs. (45). In the case where the drag coefficient, 
C d , is not being estimated and the effect of J 2 is neglected, the non-zero elements of F 0 are as 
follows 10 . Let 


_ Vp7(l + e-cos9 y 
6 [a(l-e 2 )f 2 


U =-— f 6 (X) (50.a) 

2a 


^a6,2 — 


3e 2cos0 

— -j 

1-e l + e-cos0 


f 6 (X) (50.b) 


f 2e ' sin9 f 

* 6 ' 6 1 + e • cos 0 * 6 ) 


In the more general case where C d and the influence of J 2 are considered, the development of 


F 0 is more lengthy and is found in Refs. 10 and 11. 
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Attitude Error-State Dynamics 

As mentioned earlier, the development of F a , the dynamics model that describes the 

evolution of the attitude and gyro bias errors between measurement updates, is done entirely 
differently than that of F 0 . From Eq. (37) 


dq T =[V,l] 

(51) 

Generalization of the relation expressed in Eq. (38) yields 


q = q ® dq 

from which we obtain 

(52) 

q = q®dq + q®dq 

(53) 


Noting that q" 1 ® q = u , where u is the unity quaternion, we get from the last equation 

dq = q~‘ ® q - q _1 ® q ® dq (54) 

We can replace q by yq® to, where to is to presented in quaternion form [see Ref 16], and 

similarly q by yq ® to m , where to m is the measured angular rate, to m , presented also in 
quaternion form. The last equation then becomes 

dq = q~‘ ® jq ® to m - q' 1 ® |q ® oo ® dq (55) 

Using Eq. (52) as well as the relation to m = to + dto (where dto is the angular rate measurement 
error, which, in our case, is the gyro drift rate) Eq. (55) can be written as 

dq = q“‘ ® yq ® dq ® (to + dto) - q' 1 ® yq ® w ® dq (56) 

Eq. (56) can be written as 

dq =y[dq®to-co®dq] + ydq®dco (57) 

Using the rules of quaternion multiplication [see e.g. Ref. 11, p. 759], the real part of the 
quaternion [dq ® to - to ® dq] is identically zero; that is 
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[dq 0 co - to 0 dq] r =0 


(58. a) 


where r denotes the real part. Similarly we find that 

[dq 0 dco] r = -da> ■ ya = 0 (58.b) 

where dto is a vector. For the imaginary part of the quaternion [dq 0 to - to 0 dq] we find that 

[dq <S> to- co 0 dq] Im = 2<ox ja (59. a) 

where Im denotes the imaginary part, and to is a vector. For the imaginary part of the 

quaternion [dq 0 dto] we find that 

[dq 0 dto] [m = dto + dto x « dto (59. b) 

Using Eqs. (58) and (59), Eq. (57) becomes 

a = to x a + dto (60. a) 

The term dto is the error introduced by the gyros, which we assume it consists of two parts: a 
very low frequency signal, and a wide band signal. We modeled the former as bias and the latter 
as white noise; thus 

dto = b + n (60.b) 

where b is the bias state whose model was given in Eq. (47.b). As mentioned earlier, b was 
augmented with the original state vector (see Eq. 2). In the error state vector, x, defined in Eq. 
(10), the error in estimating b is denoted 5b where obviously 

5b = 0 (60.c) 

The resulting matrix, F a , describes the propagation of a as well as the estimate of the 
gyro constant drift rate. It is a 6x6 matrix. In view of Eq. (60) the non-zero elements of F a are 
as follows. 

f l4 , =1 (61. a) f a4 .5=“z (6 lb) fa4,6=-®y ( 61C ) 
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f, s . ; - 1 

(61.d) 

U* = -to. 

(61. e) 

f.5. 6 =to x (61. g) 

= 1 

(61. h) 

f j6 .4 = 

(61. i) 

(61 j) 


HI. FILTER TESTING 


The observability of the unified filter was examined first using the following analysis. Let 
us denote the transition matrix, which corresponds to F of Eq. (48), by A k (X). This matrix 
transforms x k _, , the error in the state estimate at time t k _, , to x k the error in the state estimate at 
time t k . If at a certain time point t m , the initial error in the state estimate, denoted by x 0 , can be 
computed then for our purposes the system is observable. This is so because in the EKF, x is 
defined as the difference between the system state vector, Xand its known estimate, X. 
Adopting the common approach to the proof of complete observability of a discrete linear 
system, we express the first m measurements as follows 

r j i 


yj= H i 




x 0 j = 0, 1, m-1 


(62) 


|_i=o 

where A 0 (X) = I . Form the matrix equation 


’ y 0 " 


H 0 

y« 


H,Aj(X) 

y 2 

= 

H 2 A 2 (X)A t (X) 

_y m -i_ 


_H b . 1 A w _ i (X)...A 1 (X)_ 


(63) 


If there are n independent rows in the right hand side matrix (the observability matrix) in Eq. 
(63) then, of course, x 0 can be computed, hence the state is observable. Testing the filter with 
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CGRO data, it was found that after 5 measurements there were 13 independent rows in the 
observability matrix, and because the size of the state vector was 13, the system was observable. 

After finding that the filter operated on an observable CGRO system, the EKF algorithm 
described in the preceding was tested on data obtained from four satellites. The four satellites 
were, CGRO, RXTE, ERBS, and TOMS. The launch dates of these satellites and the start time of 
the data used in the testing is presented in Table I. The satellites varied in altitude, 


Table I. Satellite Launch Times 


Satellite 

Launch 

Date 

Start of Data 

RXTE 

Dec. 30, 1995 

1996, Sept. 8, 19h53min. 11.769sec 

ERBS 

Oct. 5, 1984 

1984, Dec. 25, 19h 22min. 12.000sec 

CGRO 

April 5, 1991 

1993, May 3, 14h 12min. 36.841 sec 

TOMS 

July 2, 1996 

1997, March 15, 19h 5min. 14.5 Msec 


inclination, resolution of the telemetry of the magnetometer data (quantization error), and 
frequency at which the magnetometer telemetry was received. Table II summarizes these 


Table II. Satellite Information 


Satellite 

h 

(km) 

i (deg) 

Error 

(mG) 

■ 

ISI 

RXTE 

580 

23.0 

0.3 

2 

ERBS 

614 

57 

6.4 

16 

CGRO 

340 

28.5 

0.3 

3-4 

TOMS 

483 

97 

4 

33-34 


variations, where h is the average altitude (the orbits were nearly circular), i is the inclination, 
‘Error’ is the quantization error of the telemetry, and AT is the time between magnetometer 
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measurements. In addition, the accuracy of the gyroscopes varied. (Note that the gyro data was 
needed for computing the attitude dynamics model of Eqs. 47. a and 61.) The CGRO and RXTE 
gyros were considerably more accurate than those of TOMS and ERBS. (CGRO used the 
Teledyne DRIRU-II gyro whose drift rate stability was 0.003 degrees per hour over six hours. 
RXTE used the Kearfott SKIRU-DII which was equivalent to the DRIRU-II gyro. TOMS and 
ERBS used less accurate gyros whose drift rate stability was about an order of magnitude larger). 
The initial estimate of each satellite position, velocity and attitude determines the initial orbital 
parameters and the initial quaternion, which constitute the initial state estimate. 


Table III. Initial RSS Errors 


Satellite 

Position 

(km) 

Velocity 

(km/sec) 

Attitude 

(deg) 

RXTE 

2617 

2.6 

15.7 

ERBS 

1000 

1.1 

11.4 

CGRO 

1098 

1.1 

12.9 

TOMS 

2389 

2.6 

14 


The initial estimation errors were determined by comparing the filter position, velocity, and 
attitude with estimates of position, velocity, and attitude computed by NASA GSFC flight 
dynamics personnel. The initial errors, which correspond to the initial estimate chosen for each 
satellite are given in Table III. The errors displayed are the RSS errors. 
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Fig. 2: The Evolution of the RXTE RSS Position Estimation Error. 
(From Initial Error of 2617 km). 



Fig. 3: The Evolution of the RXTE RSS Attitude Estimation Error. 
(From Initial Error of 15.7 deg). 


The figures presented here show the evolution of the position and attitude estimation errors for 
each of the four satellites. Figures 2 and 3 show, respectively, the RXTE RSS position and 


attitude errors, Figs. 4 and 5 show ERBS results, Figs. 6 and 7 show the CGRO results, and 
finally, TOMS results are presented in Figs. 8 and 9. The velocity errors are not shown, but for 
each satellite they have a shape similar to the position error. 




12 3 4 5 6 7 8 

time (orbit*) 


Fig. 5: The Evolution of the ERBS RSS Attitude Estimation Error. 
(From Initial Error of 1 1.4 deg). 
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Fig. 6: The Evolution of the CGRO RSS Position Estimation Error. 
(From Initial Error of 1098 km). 



Fig. 7: The Evolution of the CGRO RSS Attitude Estimation Error. 
(From Initial Error of 12.9 deg). 
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Fig. 8: The Evolution of the TOMS RSS Position Estimation Error. 
(From Initial Error of 2389 km). 



Fig. 9: The Evolution of the TOMS RSS Attitude Estimation Error. 
(From Initial Error of 14 deg). 
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Table IV, summarizes the average RSS errors for position, velocity, and attitude for each 
satellite (excluding transients). Note that the RXTE transients last for 15 orbits therefore the 
RXTE RSS errors are averaged only over the last 3 orbits. For ERBS, CGRO, and TOMS 

Table IV : Average, RSS Errors (excluding transient effects) 


Satellite 

Position 

Velocity 

Attitude 


(km) 

(km/sec) 

(deg) 

RXTE* 

15 

0.015 

1.0 

ERBS 

25 

0.03 

1.4 

CGRO 

20 

0.02 

0.2 

TOMS 

20 

0.025 

0.75 


’Averages are for the last 3 orbits. 

the final position errors are comparable, with averages of 20 to 30 km. All three also converge 
quickly, including TOMS which started with a 2389 km position error. The position estimates 
from all three appear to still be converging. However, in order to assure convergence for all four 
satellites, data of one more earth revolution inside the orbit was necessary. That is, data of some 
additional 16 orbits was required. Unfortunately the data spans that were available for these 
satellites were limited. However, the resemblance of the results to results obtained via 
simulations where, of course, data was not limited, confirm the convergence of the results with 

real data. 

CGRO has the lowest attitude error. The final average RSS attitude error is 0.2 deg. The 
CGRO data was extensively calibrated prior to its use, the magnetometer quantization error was 
small and the measurement frequency, particularly of the gyro data was high (0.256 sec between 
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gyro readings). Surprisingly, the TOMS attitude errors were also relatively small, with a final 
average of 0.75 degrees. This was a result of the high inclination for TOMS, which resulted in 
good observability, despite the large quantization error, the low measurement frequency, and the 
un-calibrated magnetometer data. ERBS had the highest attitude error, with a final average of 1.4 
deg. This was because ERBS had the largest quantization error and also a low measurement 
frequency. 

RXTE took considerably longer to converge, and both the attitude and position error results 
show a divergence in the middle of the span. It is noted that RXTE underwent three attitude 
maneuvers during this span, which apparently led to a deviation in the drag coefficient, which 
then led to the increased errors. Once the drag coefficient converged again, the position results 
improved and for the last 3 orbits were very good with an average RSS position error of 
approximately 15 km. The final average RSS attitude errors were approximately 1 degree. The 
RXTE data was not calibrated at all and contains unknown and uncompensated disturbances, 
which contributed to the larger final attitude error. 

These results compare favorably with results obtained in Ref. 17 where only magnetometers 
were used. In particular the present results fit the accuracy predicted in that reference. They also 
fit results reported in Ref. 18, which were obtained with simulated data and are better than the 
results presented in Reference 19. 

The gyros that were used on board these spacecraft were well calibrated, therefore their bias 
was negligible. Consequently the gyro bias states were deactivated; however, in order to examine 
the ability of the filter to estimate bias, an artificial bias of 1.0 mrad/sec was added to the CGRO 
gyro output and the filter was re-run with active bias states. The bias states were estimated 
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without corrupting the estimation ot the other states. Figure 10 presents the time history ot the 
estimated bias versus its correct value in the three body axes. 



Tne (Irian) 


Fig 10: The CGRO Artificial Bias Estimates as a Function of Time (Orbits) 

While successful in estimating the larger part of the gyro bias, the filter was unable to estimate 
very low gyro bias values, apparently because of the large measurement noise. Moreover, when 
the bias states were activated and the bias to be estimated was negligible, the orbit and attitude 
estimates degraded considerably. On the other hand, in the presence of large gyro bias, the 
inclusion of the bias states in the filter was imperative to obtain good orbit and attitude estimates. 
A byproduct of this inclusion was the estimation of the major part of the gyro bias. The 
degradation in performance of the filter, when the bias states of the estimator were activated 
despite the lack of actual bias, has been explained before by the Information Dilution Theorem' 0 . 
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IV. CONCLUSIONS 


This paper presented a new unified EKF algorithm tor estimating spacecraft orbit and attitude 
based on measurements of the magnetic field and the angular velocity of the spacecraft. To test 
this algorithm, it was first established that the system, which the EKF had to estimate, was 
indeed observable for at least the CGRO satellite. Next, magnetometer data and gyro data from 
four satellites were processed by the EKF. All four converged to final averages of 15 - 30 km in 
position, 0.015 - 0.03 km/sec in velocity, and 0.2 - 1.5 degrees in attitude (all RSS). Additional 
results indicated that the EKF could converge from extremely large initial position and velocity 
errors; CGRO and TOMS, for example, overcame initial position errors exceeding 5000 km. 

The data from each of the satellites differed in quantization error and measurement frequency, 
and the satellites were at 4 different inclinations. The higher inclinations for ERBS and TOMS 
resulted in quick convergence, and for TOMS gave good final results despite inaccuracies in the 
data and a large time between measurements. The coarse quantization of the ERBS telemetry and 
the infrequency of the measurements resulted in the lowest accuracy. On the other hand, the 
extensive calibration of the CGRO data, the small quantization errors, and the high data 
frequency contributed to the quick convergence and the higher accuracy achieved with CGRO. It 
was observed that RXTE required more time to converge, with good final results in position. The 
RXTE magnetometer data contained disturbances that were not calibrated and hence contributed 
to the errors in attitude. 

Although the filter was designed to also estimate gyro bias, in our experiments we deactivated 
the bias states because the available gyro data had been calibrated (which removed the gyro bias). 
In order to test the filter performance when gyro bias is present, artificial gyro bias was added to 
the CGRO data and the filter bias states were activated. It was found that the filter estimated the 
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gyro bias without degrading the quality of the attitude and orbit estimation. The filter 
successfully estimated the larger part of the gyro bias, but was unable to estimate very low gyro 
bias values, probably because of the large measurement noise. Moreover, when trying to estimate 
negligible gyro drifts, the orbit and attitude estimates degraded considerably. However, in the 
presence of large gyro bias, it was necessary to activate the bias states in the filter in order to 
obtain good orbit and attitude estimates. A byproduct of this inclusion was the estimation of the 
major part of the gyro bias. 

In summary it was shown that magnetometer and gyro measurements on board low earth- 
orbiting satellites are sufficient for coarse orbit as well as attitude determination when processed 
by the single EKF, which was developed in this work. 
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